subroutine exe_event(lunit,cunit,eid)
  use param
  use event_param
  use run_param,only:n_exe_event
  use mat_param,only:dist_r0
  use mat_param,only:half_a0
  use label_data_param,only:obj_type
  use table_param,only:event_table
  use table_param,only:sim_obj
  use geom_data_param,only:geom_obj
  use react_data_param,only:react_type
  use output_param,only:print_reaction_rate
  use group_param,only:react_group
  use sink_data_param,only:sink_type

  implicit none

  integer,intent(in) :: lunit,cunit,eid

  integer :: i
  integer :: sid,lid,bid,hid,did,xid,gid,pid,rid,tid
  integer :: sid2,lid2
  integer :: event_type,drank,dir,np,rrank
  real(dbl) :: rn,rn2,d1,d2,d3,dr
!  real(dbl) :: px,py,pz,px0,py0,pz0
  real(dbl) :: pos(3),pos0(3)
  logical :: removed

  integer :: pos1,pos2
  character(len=len_linebuf) :: str

  event_type=event_table(eid)%event_type

#ifdef DEBUG
  n_exe_event=n_exe_event+1
  call log_msg1(cunit,"")
  pos1=0
  pos2=pos1+18
  write(str(pos1+1:pos2),'("exe_event (",i6,")")') n_exe_event
  pos1=pos2
  call log_msg1_indent(cunit,1,str(1:pos2))
  call print_event(cunit,eid,1)
#endif

  select case(event_type)
  case(diff_event)

     ! Obtain the indices
     sid=event_table(eid)%obj1
     lid=sim_obj(sid)%object_type

     ! Current position
     do i=1,3
        pos(i)=sim_obj(sid)%pos(i)
     end do

     ! Determine new position
     select case(obj_type(lid)%motion_type)
     case(0)
        ! pass
     case(1)
        ! 1-D motion
        call random_number(rn)
        if (rn.le.0.5) then
           do i=1,3
              pos(i)=pos(i)+sim_obj(sid)%dr(i)
           end do
        else
           do i=1,3
              pos(i)=pos(i)-sim_obj(sid)%dr(i)
           end do
        end if
     case(2)
        ! not ready
        call abnormal_stop("exe_event",&
             "2-D motion is not ready.")
     case(3)
        ! Now it ignores the alignment of SIA.
        ! The consideration of direction requires 
        ! the energy for rotation.
        ! This wil be implemented later.
        call random_number(rn)
        dir=ceiling(rn*8.0)
        select case(dir)
        case(1)
           pos(1)=pos(1)+half_a0
           pos(2)=pos(2)+half_a0
           pos(3)=pos(3)+half_a0
        case(2)
           pos(1)=pos(1)-half_a0
           pos(2)=pos(2)+half_a0
           pos(3)=pos(3)+half_a0
        case(3)
           pos(1)=pos(1)-half_a0
           pos(2)=pos(2)-half_a0
           pos(3)=pos(3)+half_a0
        case(4)
           pos(1)=pos(1)+half_a0
           pos(2)=pos(2)-half_a0
           pos(3)=pos(3)+half_a0
        case(5)
           pos(1)=pos(1)+half_a0
           pos(2)=pos(2)+half_a0
           pos(3)=pos(3)-half_a0
        case(6)
           pos(1)=pos(1)-half_a0
           pos(2)=pos(2)+half_a0
           pos(3)=pos(3)-half_a0
        case(7)
           pos(1)=pos(1)-half_a0
           pos(2)=pos(2)-half_a0
           pos(3)=pos(3)-half_a0
        case(8)
           pos(1)=pos(1)+half_a0
           pos(2)=pos(2)-half_a0
           pos(3)=pos(3)-half_a0
        end select
     case(4)
        ! Ignore BCC lattice structure
        call random_number(rn)
        rn=rn*2.0*pi
        call random_number(rn2)
        rn2=rn2*2.0*pi
        d1=cos(rn2)
        pos(1)=pos(1)+dist_r0*cos(rn)*d1
        pos(2)=pos(2)+dist_r0*sin(rn)*d1
        pos(3)=pos(3)+dist_r0*sin(rn2)
     case default
        call abnormal_stop("exe_event",&
             "Unexpected motion type is detected.")
     end select

     call boundary_check(pos(1),pos(2),pos(3),removed)

     if (removed) then
        call remove_object(sid,lid)

        ! Measure the reaction rate
        if (print_reaction_rate) then
           do i=1,obj_type(lid)%n_free_react_group
              tid=obj_type(lid)%free_react_group(i)
              react_group(tid)%a_count=&
                   react_group(tid)%a_count+1
           end do
        end if
     else

        call move_object(lid,sid,pos(1),pos(2),pos(3))

     end if

  case(diss_event)

     ! Obtain the indices
     sid=event_table(eid)%obj1
     drank=event_table(eid)%sqnc
     lid=sim_obj(sid)%object_type

     ! Current position
     do i=1,3
        pos0(i)=sim_obj(sid)%pos(i)
     end do

     d1=sim_obj(sid)%vec(1)
     d2=sim_obj(sid)%vec(2)
     d3=sim_obj(sid)%vec(3)

     ! Remove the object
     call remove_object(sid,lid)

     ! Number of products resulting from the dissociation.
     np=obj_type(lid)%n_diss(drank)

     ! Measure the reaction rate
     if (print_reaction_rate) then
        do i=1,obj_type(lid)%n_diss_react_group
           tid=obj_type(lid)%diss_react_group(i)
           react_group(tid)%a_count=&
                react_group(tid)%a_count+1
        end do
     end if

     ! Dissociation type
     if (np.eq.0) go to 100

     ! Get heaviest type
     hid=1
     do i=2,np
        bid=obj_type(lid)%diss_obj(drank)%diss_type(hid)
        did=obj_type(lid)%diss_obj(drank)%diss_type(i)
        if (obj_type(bid)%diff_const_rate.gt.&
             obj_type(did)%diff_const_rate) hid=i
     end do

     ! Determine new type
     did=obj_type(lid)%diss_obj(drank)%diss_type(hid)

     ! Create a new object by diff_barrier at original position
     dir=0
     call create_object(did,pos0(1),pos0(2),pos0(3),&
          dir,d1,d2,d3)

     do i=1,np

        ! the heaviest object
        if (i.eq.hid) cycle

        ! the dissociated object
        did=obj_type(lid)%diss_obj(drank)%diss_type(i)
        dr=obj_type(lid)%diss_rd(drank)

        ! determine the new position
        call set_new_pos(did,dr,pos0(1),pos0(2),pos0(3),&
             d1,d2,d3,pos(1),pos(2),pos(3))

        ! create object
        dir=0
        call create_object(did,pos(1),pos(2),pos(3),&
             dir,d1,d2,d3)

     end do

  case(rctn_event)

     ! Obtain the indices for two objects
     sid=event_table(eid)%obj1
     sid2=event_table(eid)%obj2
     lid=sim_obj(sid)%object_type
     lid2=sim_obj(sid2)%object_type
     rid=event_table(eid)%react_index
     rrank=event_table(eid)%sqnc

     ! Determine the heaviest object
     if (obj_type(lid)%diff_const_rate.lt.&
          obj_type(lid2)%diff_const_rate) then
        d1=sim_obj(sid)%vec(1)
        d2=sim_obj(sid)%vec(2)
        d3=sim_obj(sid)%vec(3)
     else if (obj_type(lid)%diff_const_rate.gt.&
          obj_type(lid2)%diff_const_rate) then
        d1=sim_obj(sid2)%vec(1)
        d2=sim_obj(sid2)%vec(2)
        d3=sim_obj(sid2)%vec(3)
     else
        call random_number(rn)
        if (rn.le.0.5) then
           d1=sim_obj(sid)%vec(1)
           d2=sim_obj(sid)%vec(2)
           d3=sim_obj(sid)%vec(3)
        else
           d1=sim_obj(sid2)%vec(1)
           d2=sim_obj(sid2)%vec(2)
           d3=sim_obj(sid2)%vec(3)
        end if
     end if

     ! Currrent position
     do i=1,3
        pos0(i)=sim_obj(sid)%pos(i)
     end do

     ! remove the one object
     call remove_object(sid,lid)

     ! Current position
     do i=1,3
        pos(i)=sim_obj(sid)%pos(i)
     end do

     ! remove the other one
     call remove_object(sid2,lid2)

     ! Measure the reaction rate
     if (print_reaction_rate) then
        do i=1,react_type(rid)%n_react_group
           tid=react_type(rid)%react_group(i)
           react_group(tid)%a_count=&
                react_group(tid)%a_count+1
        end do
     end if

     ! Number of products
     np=react_type(rid)%n_new(rrank)

     ! when np=0,
     if (np.eq.0) go to 100

     ! New position
     do i=1,3
        pos0(i)=0.5*(pos(i)+pos0(i))
     end do

     ! Get the heaviest one among the products
     hid=1
     do i=2,np
        lid=react_type(rid)%new(rrank)%prod_type(hid)
        lid2=react_type(rid)%new(rrank)%prod_type(i)
        if (obj_type(lid)%diff_const_rate.gt.&
             obj_type(lid2)%diff_const_rate) hid=i
     end do

     ! The type of the heaviest product
     lid=react_type(rid)%new(rrank)%prod_type(hid)

     ! Create the object
     dir=0
     call create_object(lid,pos0(1),pos0(2),pos0(3),&
          dir,d1,d2,d3)

     ! Create the rest objects
     do i=1,np
        if (i.eq.hid) cycle
        ! Determine the new object
        lid=react_type(rid)%new(rrank)%prod_type(i)

        ! distance from the heaviest (capture radius+0.5)
        rn2=react_type(rid)%react_rd(rrank)+0.0

        ! New position
        call set_new_pos(lid,rn2,pos0(1),pos0(2),pos0(3),&
             d1,d2,d3,pos(1),pos(2),pos(3))

        ! Create the object
        dir=0
        call create_object(lid,pos(1),pos(2),pos(3),&
             dir,d1,d2,d3)

     end do

  case(extn_event)

     xid=event_table(eid)%obj1
     pid=event_table(eid)%profile_table_index

     call exe_external_event(lunit,cunit,xid,pid)

  case(sink_event)

     ! Obtain the indices
     sid=event_table(eid)%obj1
     lid=sim_obj(sid)%object_type
     rid=event_table(eid)%react_index

     ! Remove the object
     call remove_object(sid,lid)

     ! Measure the reaction rate
     if (print_reaction_rate) then
        do i=1,sink_type(rid)%n_react_group
           tid=sink_type(rid)%react_group(i)
           react_group(tid)%a_count=&
                react_group(tid)%a_count+1
        end do
     end if

  case(srce_event)

     ! Determine the new type
     gid=event_table(eid)%obj1
     rrank=event_table(eid)%obj_seq1

     np=geom_obj(gid)%diss_n(rrank)

     do i=1,np
        ! object type
        lid=geom_obj(gid)%diss_obj(rrank)%diss_type(i)

        ! Sourcejohn5!833 position
        call set_source_pos(gid,rrank,lid,pos(1),pos(2),pos(3))

        call boundary_check(pos(1),pos(2),pos(3),removed)

        if (removed) cycle

        dir=1
        call create_object(lid,pos(1),pos(2),pos(3),&
             dir,d1,d2,d3)
     end do

     ! Measure the reaction rate
     if (print_reaction_rate) then
        do i=1,geom_obj(gid)%n_diss_react_group
           tid=geom_obj(gid)%diss_react_group(i)
           react_group(tid)%a_count=&
                react_group(tid)%a_count+1
        end do
     end if

  case default
     call abnormal_stop("exe_event",&
          "Unexpected event type is detected.")
  end select

100 continue

#ifdef DEBUG
  call print_sim_obj(cunit)
  call print_event_profile(cunit)
  call print_event_table(cunit)
  call print_del_event_table(cunit)
  call print_spont_event_table(cunit)
#endif

  return

end subroutine exe_event
